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ABSTRACT 

This paper presents an algorithm for multiobjective optimization that blends together a 
number of heuristics. A population of agents combines heuristics that aim at exploring 
the search space both globally and in a neighborhood of each agent. These heuristics 
are complemented with a combination of a local and global archive. The novel agent- 
based algorithm is tested at first on a set of standard problems and then on three specific 
problems in space trajectory design. Its performance is compared against a number of 
state-of-the-art multiobjective optimisation algorithms that use the Pareto dominance 
as selection criterion: NSGA-II, PAES, MOPSO, MTS. The results demonstrate that 
the agent-based search can identify parts of the Pareto set that the other algorithms 
were not able to capture. Furthermore, convergence is statistically better although the 
variance of the results is in some cases higher. 
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1. INTRODUCTION 

The design of a space mission steps through different phases of increasing complexity. 
In the first phase, a trade-off analysis of several options is required. The trade-off 
analysis compares and contrasts design solutions according to different criteria and 
aims at selecting one or two options that satisfy mission requirements. In mathematical 
terms, the problem can be formulated as a multiobjective optimization problem. 
As part of the trade-off analysis, multiple transfer trajectories to the destination need to 
be designed. Each transfer should be optimal with respect to a number of criteria. The 
solution of the associated multiobjective optimization problem, has been addressed, by 
many authors, with evolutionary techniques. Coverstone et al. HI proposed the use 
of multiobjective genetic algorithms for the optimal design of low-thrust trajectories. 
Dachwald et al. proposed the combination of a neurocontroller and of a multiobjective 
evolutionary algorithm for the design of low-thrust trajectories [2]. In 2005 a study 
by Lee et al. |3 1 proposed the use of a Lyapunov controller with a multiobjective evo- 
lutionary algorithm for the design of low-thrust spirals. More recently, Schiitze et al. 
proposed some innovative techniques to solve multiobjective optimization problems 
for multi-gravity low-thrust trajectories. Two of the interesting aspects of the work 
of Schiitze et al. are the archiving of e- and A-approximated solutions, to the known 
best Pareto front [j?], and the deterministic pre-pruning of the search space 15]. In 
2009, Delnitz et al. ||6l proposed the use of multiobjective subdivision techniques for 
the design of low-thrust transfers to the halo orbits around the L2 libration point in the 
Earth-Moon system. Minisci et al. presented an interesting comparison between an 
EDA-based algorithm, called MOPED, and NSGA-II on some constrained and uncon- 
strained multi-impulse orbital transfer problems flT] . 

In this paper, a hybrid population-based approach that blends a number of heuristics is 
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proposed. In particular, the search for Pareto optimal solutions is carried out globally 
by a population of agents implementing classical social heuristics and more locally by 
a subpopulation implementing a number of individualistic actions. The reconstruction 
of the set of Pareto optimal solutions is handled through two archives: a local and a 
global one. 

The individualistic actions presented in this paper are devised to allow each agent to 
independently converge to the Pareto optimal set. Thus creating its own partial repre- 
sentation of the Pareto front. Therefore, they can be regarded as memetic mechanisms 
associated to a single individual. It will be shown that individualistic actions signifi- 
cantly improve the performance of the algorithm. 

The algorithm proposed in this paper is an extension of the Multi- Agent Collaborative 
Search (MACS), initially proposed in |I8]|9], to the solution of multiobjective optimisa- 
tion problems. Such an extension required the modification of the selection criterion, 
for both global and local moves, to handle Pareto dominance and the inclusion of new 
heuristics to allow the agents to move toward and along the Pareto front. As part of 
these new heuristics, this papers introduces a dual archiving mechanism for the man- 
agement of locally and globally Pareto optimal solutions and an attraction mechanisms 
that improves the convergence of the population. 

The new algorithm is here applied to a set of known standard test cases and to three 
space mission design problems. The space mission design cases in this paper consider 
spacecraft equipped with a chemical engine and performing a multi-impulse transfer. 
Although these cases are different from some of the above-mentioned examples, that 
consider a low-thrust propulsion system, nonetheless the size and complexity of the 
search space is comparable. Furthermore, it provides a first test benchmark for multi- 
impulsive problems that have been extensively studied in the single objective case but 
for which only few comparative studies exist in the multiobjective case |7|. 
The paper is organised as follows: section two contains the general formulation of the 
problem, the third section starts with a general introduction to the multi-agent collab- 
orative search algorithm and heuristics before going into some of the implementation 
details. Section four contains a set of comparative tests that demonstrates the effec- 
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tiveness of the heuristics implemented in MACS. The section briefly introduces the 
algorithms against which MACS is compared and the two test benchmarks that are 
used in the numerical experiments. It then defines the performance metrics and ends 
with the results of the comparison. 

2. PROBLEM FORMULATION 

A general problem in multiobjective optimization is to find the feasible set of solutions 
that satisfies the following problem: 

minf(x) (1) 

where D is a hyperrectangle defined as D ^ {xj \ Xj G [6^ 6"] ^ K, J = 1, 
and f is the vector function: 

f:D^R"\ f(x) = [/i(x),/2(x),...,/™(x)]^ (2) 

The optimality of a particular solution is defined through the concept of dominance: 
with reference to problem ([T), a vector y G Z? is dominated by a vector x G if 
/j (x) < fj{y) for all j — 1, ...,m. The relation x ^ y states that x dominates y. 
Starting from the concept of dominance, it is possible to associate, to each solution in 
a set, the scalar dominance index: 

Id{xj) = \{i\iAjeNpAx,-< Xj}\ (3) 

where the symbol |.| is used to denote the cardinality of a set and Np is the set of the 
indices of all the solutions. All non-dominated and feasible solutions form the set: 

X = {xeD\ Id{x) - 0} (4) 

Therefore, the solution of problem ([l) translates into finding the elements of X. If X 
is made of a collection of compact sets of finite measure in R", then once an element 
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of X is identified it makes sense to explore its neighborhood to look for other elements 
of X. On the other hand, the set of non dominated solutions can be disconnected and 
its elements can form islands in D. Hence, restarting the search process in unexplored 
regions of D can increase the collection of elements of X. 

The set X is the Pareto set and the corresponding image in criteria space is the Pareto 
front. It is clear that in D there can be more than one Xi containing solutions that are 
locally non-dominated, or locally Pareto optimal. The interest is, however, to find the 
set Xg that contains globally non-dominated, or globally Pareto optimal, solutions. 

3. MULTIAGENT COLLABORATIVE SEARCH 

The key motivation behind the development of multi-agent collaborative search was to 
combine local and global search in a coordinated way such that local convergence is 
improved while retaining global exploration 1(9] . This combination of local and global 
search is achieved by endowing a set of agents with a repertoire of actions producing 
either the sampling of the whole search space or the exploration of a neighborhood 
of each agent. More precisely, in the following, global exploration moves will be 
called collaborative actions while local moves will be called individualistic actions. 
Note that not all the individualistic actions, described in this paper, aim at exploring 
a neighborhood of each agent, though. The algorithm presented in this paper is a 
modification of MACS to tackle multiobjective optimization problems. In this section, 
the key heuristics underneath MACS will be described together with their modification 
to handle problem ([T]i and reconstruct X. 

3.1. GENERAL ALGORITHM DESCRIPTION 

A population Pq of Upop virtual agents, one for each solution vector Xj, with i = 
1, Upop, is deployed in D. The population evolves through a number of generations. 
At every generation k, the dominance index (O of each agent Xi_fc in the population 
Pfc is evaluated. The agents with dominance index Id — Q form a set Xk of non- 
dominated solutions. Hence, problem ([T]l translates into finding a series of sets Xk 
such that Xk — i- Xg for k — > kmax with kmax possibly a finite number 
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The position of each agent in D is updated through a number of heuristics. Some 
are called collaborative actions because are derived from the interaction of at least two 
agents and involve the entire population at one time. The general collaborative heuristic 
can be expressed in the following form: 

Xfc = Xfc + ^(xfc + Ufc)ufc (5) 

where depends on the other agents in the population and S* is a selection function 
which yields if the candidate point + is not selected or 1 if it is selected (see 
Section lX2] l. In this implementation a candidate point is selected if its dominance index 
is better or equal than the one of Xfe. A first restart mechanism is then implemented to 
avoid crowding. This restart mechanism is taken from (|9] and prevents the agents from 
overlapping or getting too close. It is governed by the crowding factor Wc that defines 
the minimum acceptable normaUzed distance between agents. Note that, this heuristic 
increases the uniform sampling rate of D, when activated, thus favoring exploration. 
On the other hand, by setting Wc small the agents are more directed towards local 
convergence. 

After all the collaborative and restart actions have been implemented, the resulting 
updated population Pk is ranked according to Id and split in two subpopulations: 
and P^,. The agents in each subpopulation implement sets of, so called, individualistic 
actions to collect samples of the surrounding space and to modify their current location. 
In particular, the last Upop — fenpop agents belong to P^ and implement heuristics that 
can be expressed in a form similar to Eq. (|5) but with that depends only on x^. 
The remaining fenpop agents belong to Pj. and implement a mix of actions that aim 
at either improving their location or exploring the neighborhood 7Vp(xi fe), with i = 
1, fenpop- Np{xi,k) is a hyperectangle centered in x^^fc. The intersection A^p(xi^fc)n 
D represents the local region around agent Xi^k that one wants to explore. The size of 
Np{xi,h) is specified by the value p{xi^k)' the i* edge of Np has length 2p{xi,k) max{6" — 

The agents in P^ generate a number of perturbed solutions for each fc, with s = 
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1, Smax- These solutions are collected in a local archive Ai and a dominance index 
is computed for all the elements in Ai. If at least one element y-g G Ai has 7^ = 
then Xi fe *i— Ys- If multiple elements of Ai have Id = 0, then the one with the 
largest variation, with respect to Xi_fc, in criteria space is taken. Figure [T(a)] shows 
three agents (circles) with a set of locally generated samples (stars) in their respective 
neighborhoods (dashed square). The arrows indicate the direction of motion of each 
agent. The figures shows also the local archive for the first agent Ai-^ and the target 
global archive Xg. 

The p value associated to an agent is updated at each iteration according to the rule 
devised in (j9l- Furthermore, a value s{xi,k) is associated to each agent Xi^k to specify 
the number of samples allocated to the exploration of Np{xi^k)- This value is updated 
at each iteration according to the rule devised in [9 1. 

The adaptation of p is introduced to allow the agents to self-adjust the neighborhood 
removing the need to set a priori the appropriate size of Np{xi,k)- The consequence 
of this adaptation is an intensification of the local search by some agents while the 
others are still exploring. In this respect, MACS works opposite to Variable Neigh- 
borhood Search heuristics, where the neighborhood is adapted to improve global ex- 
ploration, and differently than Basin Hopping heuristics in which the neighborhood is 
fixed. Similarly, the adaptation of s{xi,k) avoids setting a priori an arbitrary number of 
individualistic moves and has the effect of avoiding an excessive sampling of Np{xi^k) 
when p is small. The value of s(xi is initialized to the maximum number of allow- 
able individualistic moves s^ax- The value of s^ax is here set equal to the number 
of dimensions n. This choice is motivated by the fact that a gradient-based method 
would evaluate the function a minimum of n times to compute an approximation of 
the gradient with finite differences. Note that the set of individualistic actions allows 
the agents to independently move towards and within the set X, although no specific 
mechanism is defined as in [10|. On the other hand the mechanisms proposed in IfTOl 
could further improve the local search and will be the subject of a future investigation. 
All the elements of X^ found during one generation are stored in a global archive Ag. 
The elements in Ag are a selection of the elements collected in all the local archives. 
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Figure 1 : Illustration of the a) local moves and archive and b) global moves and 
archive. 

Figure \Uhj\ illustrates three agents performing two social actions that yield two sam- 
ples (black dots). The two samples together with the non-dominated solutions coming 
from the local archive form the global archive. The archive Ag is used to implement an 
attraction mechanism that improves the convergence of the worst agents (see Sections 
13.4.11 ). During the global archiving process a second restart mechanism that reinitial- 
izes a portion of the population (bubble restart) is implemented. Even this second 
restart mechanism is taken from B and avoids that, if p collapses to zero, the agent 
keeps on sampling a null neighborhood. 

Note that, within the MACS framework, other strategies can be assigned to the agents 
to evaluate their moves in the case of multiple objective functions, for example a de- 
composition technique IfTTII . However, in this paper we develop the algorithm based 
only on the use of the dominance index. 

3.2. COLLABORATIVE ACTIONS 

Collaborative actions define operations through which information is exchanged be- 
tween pairs of agents. Consider a pair of agents Xi and X2, with Xi -< X2. One of the 
two agents is selected at random in the worst half of the current population (from the 
point of view of the property Id), while the other is selected at random from the whole 
population. Then, three different actions are performed. Two of them are defined by 
adding to Xi a step defined as follows: 

Ufc = ar*(x2 -xi), (6) 

and corresponding to: extrapolation on the side of xi (a = — 1, i = 1), with the further 
constraint that the result must belong to the domain D (i.e., if the step Ufc leads out of 
D, its size is reduced until we get back to D); interpolation (a = 1), where a random 
point between Xi and X2 is sampled. In the latter case, the shape parameter t is defined 
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Algorithm 1 Main MACS algorithm 



1: Initialize a population Pq of Upop agents in D, k = 0, number of function evalu- 
ations Uevai = 0> maximum number of function evaluations N,., crowding factor 

Wc 

2: for alH = 1, npop do 

3: Xi.k = Xi^k + S{xi^k + Ui,k)ui^k 

4: end for 

5: Rank solutions in Pk according to Id 

6: Re-initialize crowded agents according to the single agent restart mechanism 

7: for all i = feUpop, ripop do 

8: Generate np mutated copies of x;.fe. 

9: Evaluate the dominance of each mutated copy yj, against x^^fc, with p = 

1 , . . . , Hp 

10: if 3p|yp < Xj,fe tiien 

11: p = argmaxp IIyp - Xj,fc|| 

12: x^_fe -{- yp 

13: end if 

14: end for 

15: for all i = 1, feUpop do 

16: Generate s < Smax individual actions Ug such that = Xj^/j -|- Ug 

17: if 3s|ys -< Xj^/s then 

18: s = argmaxg Hy^ - Xj,fc|| 

19: Xi^k Ys 

20: end if 

21: Store candidate elements y^ in the local archive Ai 
22: Update p{xi^k) and s(xj^fe) 

23: end for 

24: Form Pk = Pl[j Pfe" and Ag = Ag\J Ai [j Pk 

25: Compute Id of all the elements in Ag 

26: Ag = {x|x eAgA Id{x) = A ||x - \\ > Wc} 

27: Re-initiaUze crowded agents in Pk according to the second restart mechanism 

28: Compute attraction component to Ag for all Xi^k S Pfe \ Xk 

29: k = k + 1 

30: Termination Unless nevai > -^e. GoTo Step 2 



as follows: 

t = 0.75 ^(^1) ^ ^ (7) 

^max 

The rationale behind this definition is that we are favoring moves which are closer to 
the agent with a higher fitness value if the two agents have the same s value, wlule in 
the case where the agent with highest fitness value has a s value much lower than that 
of the other agent, we try to move away from it because a small s value indicates that 
improvements close to the agent are difficult to detect. 
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The third operation is the recombination operator, a single-point crossover, where, 
given the two agents: we randomly select a component j; split the two agents into 
two parts, one from component 1 to component j and the other from component j + I 
to component n; and then we combine the two parts of each of the agents in order 
to generate two new solutions. The three operations give rise to four new samples, 
denoted by Ji, y3, y4. Then, Id is computed for the set X2, y^, y2, ya, y4. The 
element with Id = becomes the new location of X2 in D. 

3.3. INDIVIDUALISTIC ACTIONS 

Once the collaborative actions have been implemented, each agent in is mutated 
a number of times: the lower the ranking the higher the number of mutations. The 
mutation mechanisms is not different from the single objective case but the selection is 
modified to use the dominance index rather than the objective values. 
A mutation is simply a random vector such that x^ fc + G D. All the mutated 
solution vectors are then compared to x^ fc, i.e. Id is computed for the set made of 
the mutated solutions and x^ fc. If at least one element of the set has Id = then 
Xi,fc <— Yp- If multiple elements of the set have Id ~ 0, then the one with the largest 
variation, with respect to x^ fc, in criteria space is taken. 

Each agent in P^, performs at most Smax of the following individuaUstic actions: in- 
ertia, differential, random with line search. The overall procedure is summarized in 
Algorithmic and each action is described in detail in the following subsections. 

3.3.1. Inertia 

If agent i has improved from generation fc — 1 to generation k, then it follows the 
direction of the improvement (possibly until it reaches the border of D), i.e., it performs 
the following step: 

= x,,fc + AA/ (8) 
where A = min{l, max{A : ys £ D}} and A/ — (x^ fc — Hi^k-i)- 
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3.3.2. Differential 

This step is inspired by Differential Evolution |,12j . It is defined as follows: let x^^ ^jX^^fe, ^ 
be three randomly selected agents; then 

Ys = ^i,k + e [xi^.k + F{xi^^k - Xi2,fe)] (9) 

with e a vector containing a random number of and 1 (the product has to be intended 
componentwise) with probability 0.8 and F = 0.8 in this implementation. For every 
component ys[j] of y-g that is outside the boundaries defining D then ys[j] — r{h'^ — 
b'j) + b'j, with r E [7(0, 1). Note that, although this action involves more than one 
agent, its outcome is only compared to the other outcomes coming from the actions 
performed by agent x^ ^ and therefore it is considered individualistic. 

3.3.3. Random with Line Search 

This move realizes a local exploration of the neighborhood A^(xi.fe). It generates a first 
random sample ys E iV(xjjj). Then if is not an improvement, it generates a second 
sample Ys+i by extrapolating on the side of the better one between and fc: 

Ys+i = Xi^fe + A [a2r*(y^ - Xi^k) + cei{y^ - Xi^k)] (10) 

with A = min{l,max{A : ys+i G D}} and where ai,a2 E { — 1,0,1}, r E 
U{0,1) and t is a shaping parameter which controls the magnitude of the displacement. 
Here we use the parameter values ai = 0, 02 = ^ 1, ^ = 1, which corresponds to 
extrapolation on the side of X; fc, and ai = 02 = 1, t = 1, which corresponds to 
extrapolation on the side of y^. 

The outcome of the extrapolation is used to construct a second order one-dimensional 
model of Id- The second order model is given by the quadratic function /;(cr) = 
aia^+a2CT+a3 where is a coordinate along the y^+i — ys direction. The coefficients 
ai, 02 and 03 are computed so that /; interpolates the values Id{ys), IdiVs+i) and 
Id{xi,k)- In particular, for = /( = IdiVs), for ct = 1 /; Id{ys+i) and for 
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<y = {^i,k - ys)/||xi,fe - Ys II /; = Id{^i,k)- Then, a new sample ys+2 is taken at the 
minimum of the second-order model along the ys+i — Ys direction. 

Algorithm 2 Individual Actions in Pj. 

1: s = 1, stop = 

2: if Xj,fe -< Xi^k-i then 

Ys =3,fe + A(xi_fe — Xj,fe_i) 

with A = niin{l, max{A : y,, G £*}}• 
3: end if 

4: if Xj,fc -< then 
s = s + 1 

= e[xi^k - (xii.fc + (x.,,,fc - x.,2,fc))] 

Vj|y«(i) ^ D, y,{j) = r{buij) - 6,(j)) + 

with J = 1, n and r G ?7(0, 1) 
5: else stop = 1 
6: end if 

7: if Xj,fe ^ then 
s = s + l 

Generate y, G Np{xi^k)- _ 
Compute y^+i = Xj^fe + Xr*{x,^k - yj 
with A = min{l, max{A : y^ € D}} 
andreUiOA). 

Compute y,+2 = Acr„i„(y^+i -ys)/ll(y5+i -yJII' 

with_(T„j„ = arg miner {aicr^ + a2(T + Id{ys)}, 

and A = min{l, max{A : yg_,_2 € D}}. 

s = s + 2 
8: else stop = 1 
9: end if 

10: Termination Unless s > Smax or stop = 1 , GoTo Step 4 



The position of Xj^fc in is then updated with the ys that has Id = and the longest 
vector difference in the criteria space with respect to Xj fe. The displaced vectors ys 
generated by the agents in are not discarded but contribute to a local archive Ai, 
one for each agent, except for the one selected to update the location of Xj^fe. In order 
to rank the ys, the following modified dominance index is used: 



id{^^,k)= {j \ fjiys) = fji^i,k)} 
{j I /j(yJ > fj{^i,k)} 



(11) 



where k is equal to one if there is at least one component of f(xjjt) = [/i, /2, fm]'^ 
which is better than the corresponding component of f(ys), and is equal to zero other- 
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wise. 

Now, if for the s*'' outcome, the dominance index in Eq. (fTTT i is not zero but is lower 
than the number of components of the objective vector, then the agent fe is only 
partially dominating the s* outcome. Among all the partially dominated outcomes 
with the same dominance index we select the one that satisfies the condition: 

min((f(x,,fc))-f(yj),e) (12) 



where e is the unit vector of dimension m,e— t^^^'^ui:-^! . All the non-dominated and 
selected partially dominated solutions form the local archive Ai . 



3.4. THE LOCAL AND GLOBAL ARCHIVES Al AND Aa 

Since the outcomes of one agent could dominate other agents or the outcomes of other 
agents, at the end of each generation, every Ai and the whole population Pk are added 
to the current global archive Ag. The global Ag contains Xk, the current best esti- 
mate of Xg. The dominance index in Eq.Q is then computed for all the elements 
in Ag ~ Ag IJ; Al [J Pfe and only the non-dominated ones with crowding distance 
||x,i_fc — x^g II > Wc are preserved (where x^i^ is an element of Ag). 

3.4.1. Attraction 

The archive Ag is used to direct the movements of those agents that are outside Xk- 
All agents, for which Id ^ at step k, are assigned the position of the elements in Ag 
and their inertia component is recomputed as: 

A/ = r(xA, - x,,fc) (13) 

More precisely, the elements in the archive are ranked according to their reciprocal dis- 
tance or crowding factor Then, every agent for which Id ^ picks the least crowded 
element x^ not already picked by any other agent. 
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3.5. STOPPING RULE 

The search is stopped when a prefixed number N^. of function evaluations is reached. 
At termination of the algorithm the whole final population is inserted into the archive 
A,. 

4. COMPARATIVE TESTS 

The proposed optimization approach was implemented in a software code, in Matlab, 
called MACS. In previous works ||8]|9l, MACS was tested on single objective optimiza- 
tion problems related to space trajectory design, showing good performances. In this 
work, MACS was tested at first on a number of standard problems, found in literature, 
and then on three typical space trajectory optimization problems. 
This paper extends the results presented in Vasile and Zuiani (2010) lfT3l by adding 
a broader suite of algorithms for multiobjective optimization to the comparison and a 
different formulation of the performance metrics. 

4. 1 . TESTED ALGORITHMS 

MACS was compared against a number of state-of-the-art algorithms for multiobjective 
optimization. For this analysis it was decided to take the basic version of the algorithms 
that is available online. Further developments of the basic algorithms have not been 
considered in this comparison and will be included in future works. Note that, all the 
algorithms selected for this comparative analysis use Pareto dominance as selection 
criterion. 

The tested algorithms are: NSGA-II OH, MOPSO |15|, PAES US and MTS 113. A 
short description of each algorithm with their basic parameters follows. 

4.LL NSGA-II 

The Non-Dominated Sorting Genetic Algorithm (NSGA-II) is a genetic algorithm 
which uses the concept of dominance class (or depth) to rank the population. A crowd- 
ing factor is then used to rank the individuals within each dominance class. Optimiza- 
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tion starts from a randomly generated initial population. The individuals in the popula- 
tion are sorted according to their level of Pareto dominance with respect to other indi- 
viduals. To be more precise, a fitness value equal to 1 is assigned to the non-dominated 
individuals. Non-dominated individuals form the first layer (or class). Those individu- 
als dominated only by members of the first layer form the second one and are assigned 
a fitness value of 2, and so on. In general, for dominated individuals, the fitness is 
given by the number of dominating layers plus 1 . A crowding factor is then assigned 
to each individual in a given class. The crowding factor is computed as the sum of 
the Euclidean distances, in criteria space, with respect to other individuals in the same 
class, divided by the interval spanned by the population along each dimension of the 
objective space. Inside each class, the individuals with the higher value of the crowding 
parameter obtain a better rank than those with a lower one. 

At every generation, binary tournament selection, recombination, and mutation oper- 
ators are used to create an offspring of the current population. The combination of 
the two is then sorted according to dominance first and then to crowding. The non- 
dominated individuals with lowest crowding factor are then used to update the popula- 
tion. 

The parameters to be set are the size of the population, the number of generations, the 
crossover and mutation probabihty, Pc and Pm, and distribution indexes for crossover 
and mutation, r/c and r/^. respectively. Three different ratios between population size 
and number of generations were considered: 0.08, 0.33 and 0.75. The values Pc and 
Pm were set to 0.9 and 0.2 respectively and kept constant for all the tests. The values, 
5, 10 and 20, were considered for r]c, while tests were run for values of rjm equal to 5, 
25 and 50. 

4.1.2. MOPSO 

MOPSO is an extension of Particle Swarm Optimization (PSO) to multiobjective prob- 
lems. Pareto dominance is introduced in the selection criteria for the candidate solu- 
tions to update the population. MOPSO features an external archive which stores all the 
non-dominated solutions and at the same time is used to guide the search process of the 
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swarm. This is done by introducing the possibility to direct the movement of a particle 
towards one of the less crowded solutions in the archive. The solution space is subdi- 
vided into hypercubes through an adaptive grid. The solutions in the external archive 
are thus reorganized in these hypercubes. The algorithm keeps track of the crowding 
level of each hypercube and promotes movements towards less crowded areas. In a 
similar manner, it also gives priority to the insertion of new non-dominated individuals 
in less crowded areas if the extemal archive has already reached its predefined maxi- 
mum size. For MOPSO three different ratios between population size and number of 
generations were tested: 0.08 0.33 0.75. It was also tested with three different numbers 
of subdivisions of the solution space: 10, 30, and 50. The inertia component in the 
motion of the particles was set to 0.4 and the weights of the social and individualistic 
components were se to 1 

4.1.3. PAES 

PAES is a (1 H- 1) Evolution Strategy with the addition of an external archive to store 
the current best approximation of the Pareto front. It adopts a population of only a 
single chromosome which, at every iteration, generates a mutated copy. The algorithm 
then preserves the non-dominated one between the parent and the candidate solution. If 
none of the two dominates the other, the algorithm then checks their dominance index 
with respect to the solutions in the archive. If also this comparison is inconclusive, 
then the algorithm selects the one which resides in the less crowded region of the 
objective space. To keep track of the crowding level of the objective space, the latter is 
subdivided in an n-dimensional grid. Every time a new solution is added to (or removed 
from) the archive, the crowding level of the corresponding grid cell is updated. PAES 
has two main parameters that need to be set, the number of subdivisions in the space 
grid and the mutation probabiUty. Values of 1,2 and 4 were used for the former and 0.6 
0.8 and 0.9 were used for the latter. 
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4.1.4. MTS 

MTS is an algorithm based on pattern search. The algorithm first generates a population 
of uniformly distributed individuals. At every iteration, a local search is performed by 
a subset of individuals. Three different search patterns are included in the local search: 
the first is a search along the direction of each decision variable with a fixed step length; 
the second is analogous but the search is limited to one fourth of all the possible search 
directions; the third one also searches along each direction but selects only solutions 
which are evenly spaced on each dimension within a predetermined upper bound and 
lower bound. At each iteration and for each individual, the three search patterns are 
tested with few function evaluations to select the one which generates the best candidate 
solutions for the current individual. The selected one is then used to perform the local 
search which will update the individual itself. The step length along each direction, 
which defines the size of the search neighbourhood for each individual, is increased if 
the local search generated a non-dominated child and is decreased otherwise. When 
the neighbourhood size reaches a predetermined minimum value, it is reset to 40% of 
the size of the global search space. The non-dominated candidate solutions are then 
used to update the best approximation of the global Pareto front. MTS was tested with 
a population size of 20, 40 and 80 individuals. 

4.2. PERFORMANCE METRICS 

Two metrics were defined to evaluate the performance of the tested multiobjective op- 
timizers: 




where Mp is the number of elements, with objective vector g, in the true global Pareto 
front and Np is the number of elements, with objective vector f , in the Pareto front that 
a given algorithm is producing. Although similar, the two metrics are measuring two 
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different things: Mgpr is the sum, over all the elements in the global Pareto front, of the 
minimum distance of all the elements in the Pareto front Np from the the z*'' element in 
the global Pareto front. Mconv, instead, is the sum, over all the elements in the Pareto 
front Np, of the minimum distance of the elements in the global Pareto front from the 
i*'' element in the Pareto front Np. 

Therefore, if Np is only a partial representation of the global Pareto front but is a 
very accurate partial representation, then metric Mspr would give a high value and 
metric M^onv a low value. If both metrics are high then the Pareto front Np is partial 
and poorly accurate. The index Mconv is similar to the mean Euclidean distance 02, 
although in Mconv the Euclidean distance is normalized with respect to the values of the 
objective functions, while Mspr is similar to the generational distance ifTSI . although 
even for Mspr the distances are normalized. 

Given n repeated runs of a given algorithm, we can define two performance indexes: 
Pconv — P{Mconv < tolconv) or the probability that the index Mconv achieves a value 
less than the threshold tolconv and pspr — P{Mspr < tolspr) or the probability that 
the index Mspr achieves a value less than the threshold tolconv 
According to the theory developed in QUI], 200 runs are sufficient to have a 95% 
confidence that the true values of Pconv and pspr are within a ±5% interval containing 
their estimated value. 

Performance index (fl4l and (flSl l tend to uniformly weigh every part of the front. This 
is not a problem for index (fl4t but if only a relatively small portion of the front is 
missed the value of performance index (flSt might be only marginally affected. For this 
reason, we slightly modified the computation of the indexes by taking only the Mp and 
Np solutions with a normalised distance in criteria space that was higher than 10^^. 
The global fronts used in the three space tests were built by taking a selection of about 
2000 equispaced, in criteria space, nondominated solutions coming from all the 200 
runs of all the algorithms. 
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4.3. PRELIMINARY TEST CASES 

For the preliminary tests, two sets of functions, taken from the literature 01411151 . were 
used and the performance of MACS was compared to the results in fl4\ and ifTSl . 
Therefore, for this set of tests, MTS was not included in the comparison. The function 
used in this section can be found in Table[T] 

Table 1 : Multiobjective test functions 

The first three functions were taken from llTSl . Test cases Deb and Scha are two 
examples of disconnected Pareto fronts, Deb2 is an example of problem presenting 
multiple local Pareto fronts, 60 in this two dimensional case. 

The last three functions are instead taken from fi^. Test case ZDT2 has a concave 
Pareto front with a moderately high-dimensional search space. Test case ZDT6 has 
a concave Pareto front but because of the irregular nature of /i there is a strong bias 
in the distribution of the solutions. Test case, ZDT4, with dimension 10, is commonly 
recognized as one of the most challenging problems since it has 21^ different local 
Pareto fronts of which only one corresponds to the global Pareto-optimal front. 
As a preliminary proof of the effectiveness of MACS, the average Euclidean distance 
of 500 uniformly spaced points on the true optimal Pareto front from the solutions 
stored in Ag by MACS was computed and compared to known results in the literature. 
MACS was run 20 times to have a sample comparable to the one used for the other 
algorithms. The global archive was limited to 200 elements to be consistent with lfT4ll . 
The value of the crowding factor Wc, the threshold ptoi and the convergence Pmin were 
kept constant to le-5 in all the cases to provide good local convergence. 
To be consistent with ifTSl . on Deb, Scha and, Deb2, MACS was run respectively for 
4000, 1200 and 3200 function evaluations. Only two agents were used for these lower 
dimensional cases, with /e = 1/2. On test cases ZDT2, ZDTA and ZDT6, MACS 
was run for a maximum of 25000 function evaluations to be consistent with [141 , with 
three agents and fe ^ 2/3 for ZDT2 and four agents and /e = 3/4 on ZDTA and 
ZDT6. 
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The results on Deb, Scha and, Deb2 can be found in Table |2] while the results on 
ZDT2, ZDTA and ZDT6 can be found in[3l 

On all the smaller dimensional cases MACS performs comparably to MOPSO and bet- 
ter than PAES. It also performs than NSGA-II on Deb and Deb2. On Scha MACS 
performs apparently worse than NSGA-II, although after inspection one can observe 
that all the elements of the global archive Ag belong to the Pareto front but not uni- 
formly distributed, hence the higher value of the Euclidean distance. On the higher 
dimensional cases, MACS performs comparably to NSGA-II on ZDT2 but better than 
all the others on ZDT4 and ZDT6. Note in particular the improved performance on 
ZDT4. 

Table 2: Comparison of the average Euclidean distances between 500 uniformly 
space points on the optimal Pareto front for various optimization algorithms: 
smaller dimension test problems. 



Table 3: Comparison of the average Euclidean distances between 500 uniformly 
space points on the optimal Pareto front for various optimization algorithms: 
larger dimension test problems. 

On the same six functions a different test was run to evaluate the performance of differ- 
ent variants of MACS. For all variants, the number of agents, /e, Wc, Ptoi and pmin was 
set as before, but instead of the mean Euclidean distance, the success rates Pconv and 
Pspr were measured for each variant. The number of function evaluations for ZDT2, 
ZDTA, ZDT6 and Deb2 is the same as before, while for Scha and Deb it was re- 
duced respectively to 600 and 1000 function evaluations given the good performance 
of MACS already for this number of function evaluations. Each run was repeated 200 
times to have good confidence in the values of Pconv and Pspr. 

Four variants were tested and compared to the full version of MACS. Variant MACS no 
local does not implement the individualistic moves and the local archive, variant MACS 
p = 1 has no adaptivity on the neighborhood, its size is kept fixed to 1, variant MACS 
p ~ 0.1 has the size of the neighborhood fixed to 0.1, variant MACS no attraction has 
the attraction mechanisms not active. 
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The result can be found in Table |4] The values of tolconv and tolgpr are respectively 
0.001 and 0.0035 for ZDT2, 0.003 and 0.005 for ZDTi, 0.001 and 0.025 for ZDT6, 
0.0012 and 0.035 for Deb, 0.0013 and 0.04 for Scha, 0.0015 and 0.0045 for Deb2. 
This thresholds were selected to highlight the differences among the various variants. 
The table shows that the adaptation mechanism is beneficial in some cases although, 
in others, fixing the value of p might be a better choice. This depends on the problem 
and a general rule is difficult to derive at present. Other adaptation mechanisms could 
further improve the performance. 

The use of individualistic actions coupled with a local archive is instead fundamental, 
so is the use of the attraction mechanism. Note, however, how the attraction mecha- 
nism penalizes the spreading on biassed problems like ZDTQ, this is expected as it 
accelerates convergence. 

Table 4; Comparison of different variants of MACS. 
4.4. APPLICATION TO SPACE TRAJECTORY DESIGN 

In this section we present the application of MACS to three space trajectory prob- 
lems: a two-impulse orbit transfer from a Low Earth Orbit (LEO) to a high-eccentricity 
Molniya-like orbit, a three-impulse transfer from a LEO to Geostationary Earth Orbit 
(GEO) and a multi-gravity assist transfer to Saturn equivalent to the transfer trajectory 
of the Cassini mission. The first two cases are taken from the work of Minisci et al. Q. 
In the two-impulse case, the spacecraft departs at time to from a circular orbit around 
the Earth (the gravity constant is fj.E = 3.986010^ km^^s^^) with radius tq = 6721 
km and at time tf is injected into an elliptical orbit with eccentricity ct = 0.667 
and semimajor axis gt = 26610 km. The transfer arc is computed as the solution 
of a Lambert's problem ll20l and the objective functions are the transfer time T — 
tf ~ to and the sum of the two norms of the velocity variations at the beginning and 
at the end of the transfer arc Avtot- The objectives are functions of the solution vector 
X = [to tff^ E D C The search space D is defined by the following intervals 
to e [0 10.8], and t/ G [0.03 10.8]. 

21 



In the three-impulse case, the spacecraft departs at time Iq from a circular orbit around 
the Earth with radius rg = 7000 km and after a transfer time T — ti + t2 is injected 
into a circular orbit with radius rf — 42000. An intermediate manoeuvre is performed 
at time to + ti and at position defined in polar coordinates by the radius ri and the 
angle 6i. The objective functions are the total transfer time T and the sum of the 
three impulses Avtot. The solution vector in this case is x = [to,ti,ri,9i,tf]'^ G 
D C M.^. The search space D is defined by the following intervals G [0 1.62], 
h e [0.03 21.54], ri G [7010 105410], 9i e [0.01 27r - 0.01], and e [0.03 21.54]. 
The Cassini case consists of 5 transfer arcs connecting a departure planet, the Earth, 
to the destination planet, Saturn, through a sequence of swing-by's with the planets: 
Venus, Venus, Earth, Jupiter Each transfer arc is computed as the solution of a Lam- 
bert's problem ||2T1 given the departure time from planet Pi and the arrival time at 
planet P^+i. The solution of the Lambert's problems yields the required incoming 
and outgoing velocities at each swing -by planet Vin and Vrout- The swing-by is mod- 
eled through a linked-conic approximation with powered maneuvers f22], i.e., the mis- 
match between the required outgoing velocity Vrout and the achievable outgoing ve- 
locity Vaout is compensated through a Av maneuver at the pericenter of the gravity 
assist hyperbola. The whole trajectory is completely defined by the departure time to 
and the transfer time for each leg Ti, with i — 1, 5. The normalized radius of the 
pericenter ^ of each swing-by hyperbola is derived a posteriori once each powered 
swing-by manoeuvre is computed. Thus, a constraint on each pericenter radius has to 
be introduced during the search for an optimal solution. In order to take into account 
this constraint, one of the objective functions is augmented with the weighted violation 
of the constraints: 

4 4 

/(x) = Ai>o + ^ Auj + Avf + ^ Wi(rp^j - rprnin,if (16) 

1=1 i=l 

for a solution vector x = [to, Ti, T2, Ta, T4, T^]^ . The objective functions are, the to- 
tal transfer time T ^Y/.^i and /(x). The minimum normaUzed pericenter radii are 
rpmin,i = 1.0496, rp,nin,2 = 1.0496, rprnin,3 = 1 .0627 and rp„im,4 = 9.3925. The 
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seai-ch space D is defined by the following intervals: to E [-1000, 0]MJD2000, Ti 6 
[30,400]d,T2 e [100,470]d,r3 e [30,400]d, r4 e [400, 2000] d, e [1000, 6000]d. 
The best known solution for the single objective minimization of /(x) is fbest — 
4.9307 km/s, with xbe^t = [-789.753, 158.2993,449.3859,54.7060,1024.5896,4552.7054]'^. 

4.4.1. Test Results 

For this second set of tests, each algorithm was run for the same number of function 
evaluations. In particular, consistent with the tests performed in the work of Minisci et 
al., Q we used 2000 function evaluations for the two-impulse case and 30000 for the 
three-impulse case. For the Cassini case, instead, the algorithms were run for 180000, 
300000 and 600000 function evaluations. 

Note that the version of all the algorithms used in this second set of tests is the one that 
is freely available online, written in c/c-n-. We tried in all cases to stick to the available 
instructions and recommendations by the author to avoid any bias in the comparison. 
The thresholds values for the two impulse cases was taken from [7 | and is tolconv — 
0.1, tolspr ~ 2.5. For the three-impulse case instead we considered tolconv — 5.0, 
tolspr = 5.0. For the Cassini case we used tolconv — 0.75, tolgpr = 5, instead. These 
values were selected after looking at the dispersion of the results over 200 runs. Lower 
values would result in a zero value of the performance indexes of all the algorithms, 
which is not very significant for a comparison. 

MACS was tuned on the three-impulse case. In particular, the crowding factor Wc, the 
threshold ptoi and the convergence pmm were kept constant to le-5, which is below 
the required local convergence accuracy, while fc and Upop were changed. A value of 
le-5 is expected to provide good local convergence and good density of the samples 
belonging to the Pareto front. Table |5]reports the value of performance indexes Pconv 
and pspr over 200 runs of MACS with different settings. The index pspr and the index 
Pconv have different, almost opposite, trends. However, it was decided to select the 
setting that provides the best convergence,i.e. Upop = 15 and fc = 1/3. This setting 
will be used for all the tests in this paper. 

On top of the complete algorithm, two variants of MACS were tested: one without 
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individualistic moves and local archive, denoted as no local in the tables, and one with 
no attraction towards the global archive Ag, denoted as no att in the tables. Only these 
two variants are tested on these cases as they displayed the most significant impact in 
the previous standard test cases and more importantly were designed specifically to 
improve performances. 

NSGA-II, PAES, MOPSO and MTS were tuned as well on the three-impulse case. In 
particular, for NSGA-II the best result was obtained for 150 individuals and can be 
found in Table |6l A similar result could be obtained for MOPSO, see Table |2l For 
MTS only the population was changed while the number of individuals performing 
local moves was kept constant to 5. The results of the tuning of MTS can be found in 
Table|9] For the tuning of PAES the results can be found in Table[8] 
All the parameters tuned in the three impulse case were kept constant except for the 
population size of NSGA-II and MOPSO. The size of the population of NSGA-II and 
MOPSO was set to 100 and 40 respectively on the two impulse case and was increased 
with the number of function evaluations in the Cassini case. In particular for NSGA- 
II the following ratios between population size and number of function evaluations 
was used: 272/180000, 353/300000, 500/600000. For MOPSO the following ratios 
between population size and number of function evaluations was used: 224/180000, 
447/300000, 665/600000. This might not be the best way to set the population size for 
these two algorithms but it is the one that provided better performance in these tests. 
Note that the size of the global archive for MACS was constrained to be lower than 
the size of the population of NSGA-II, in order to avoid any bias in the computation of 

The performance of all the algorithms on the Cassini case can be found in Table fT2] for 

a variable number of function evaluations. 

Figure 2: Three-impulse test case: a) Complete Pareto front, b) close-up of the 
Pareto fronts 

Table 5: Indexes pconv and pspr for different settings of MACS 
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Figure 3: Cassini test case: a) Complete Pareto front, b) close-up of the Pareto 
fronts 

For the three-impulse case, MACS was able to identify an extended Pareto front (see 
Fig |2(a)| and Fig. |2(b)| where all the non-dominated solutions from all the 200 runs are 
compared to the global front), compared to the results in [7 |. The gap in the Pareto 
front is probably due to a hmited spreading of the solutions in that region. Note the 
cusp due the transition between the condition in which 2-impulse solutions are optimal 
and the condition in which 3-impulse solutions are optimal. 

Table [TT] summarizes the results of all the tested algorithms on the two-impulse case. 
The average value of the performance metrics is reported with, in brackets, the associ- 
ated variance over 200 runs. The two-impulse case is quite easy and all the algorithms 
have no problems identifying the front. However, MACS displays a better convergence 
than the other algorithms while the spreading of MOPSO and MTS is superior to the 
one of MACS. 

The three-impulse case is instead more problematic (see Table [TOll. NSGA-II is not 
able to converge to the upper-left part of the front and therefore the convergence is 
and the spreading is comparable to the one of MACS. All the other algorithms perform 
poorly with a value of almost for the performance indexes. This is mainly due to the 
fact that no one can identify the vertical part of the front. 

Note that the long tail identified by NSGA-II is actually dominated by two points be- 
longing to the global front (see |2(b)| i. 

Table fT2]reports the statistics for the Cassini problem. On top of the performance of the 
three variants tested on the other two problems, the table reports also the result for 10 
agents and /e = 5. For all numbers of function evaluations MACS has better spreading 
than NSGA-II because NSGA-II converges to a local Pareto front. Nonetheless NSGA- 
II displays a more regular behavior and a better convergence for low number of function 
evaluations although it never reaches the best front. The Pareto fronts are represented 
in Fig. |3(a)| and Fig. |3(b)| (even in this case the figures represent all the non-dominated 
solutions coming from all the 200 runs). Note that, the minimum / returned by MACS 
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is the best known solution for the single objective version of this problem (see ||9]). All 
the other algorithms perform quite poorly on this case. 

Of all the variants of MACS tested on these problems the full complete one is per- 
forming the best. As expected, in all three cases, removing the individualistic moves 
severely penalizes both convergence and spreading. It is interesting to note that re- 
moving the attraction towards the global front is also comparatively bad. On the two 
impulse case it does not impact the spreading but reduces the convergence, while on 
the Cassini case, it reduces mean and variance but the success rates are zero. The ob- 
servable reason is that MACS converges slower but more uniformly to a local Pareto 
front. 

Finally, it should be noted that mean and variance seem not to capture the actual per- 
formance of the algorithms. In particular they do not capture the ability to identify the 
whole Pareto front as the success rates instead do. 

Table 6: NSGAII tuning on the 3-impulse case 

Table 7: PAES tuning on the 3-impulse case 
5. CONCLUSIONS 

In this paper we presented a hybrid evolutionary algorithm for multiobjective opti- 
mization problems. The effectiveness of the hybrid algorithm, implemented in a code 
called MACS, was demonstrated at first on a set of standard problems and then its per- 
formance was compared against NSGA-II, PAES, MOPSO and MTS on three space 
trajectory design problems. The results are encouraging as, for the same computational 
effort (measured in number of function evaluations, MACS was converging more accu- 

Table 8; MOPSO tuning on the 3-impulse case. 
Table 9: MTS tuning on the 3-impulse 
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rately than NSGA-II on the two-impulse case and managed to find a previously undis- 
covered part of the Pareto front of the three-impulse case. As a consequence, on the 
three-impulse case, MACS, has better performance metrics than the other algorithms. 
On the Cassini case NSGA-II appears to converge better to some parts of the front al- 
though MACS yielded solutions with better / and identifies once more a part of the 
front that NSGA-11 cannot attain. PAES and MTS do not perform well on the Cassini 
case, while MOPSO converges well locally but, with the settings employed in this 
study, yielded a very poor spreading. 

From the experimental tests in this paper we can argue that the following mechanisms 
seem to be particularly effective: the use of individual local actions with a local archive 
as they allow the individuals to move towards and within the Pareto set; the use of an 
attraction mechanism as it accelerates convergence. 

Finally it should be noted that all the algorithms tested in this study use the Pareto 
dominance as selection criterion. Different criteria, like the decomposition in scalar 
subproblems, can be equality implemented in MACS, without disrupting its working 
principles, and lead to different performance results. 
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Table 10: Summary of metrics Mconv and Mgpr, and associated peri'ormance 
indexes Pconv and Pspr, on the three impulse test cases. 



Table 11: Metrics Mconv and M^pr, and associated performance indexes pc 
and Pspr, on the two impulse test cases. 



Table 12: Metrics Mconv and Mgpr, and associated performance indexes Pc 
and Pspr, on the Cassini case. 
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NOMENCLATURE 





global archive 


Ai 


local archive 




semimajor axis 


D 


search space 




eccentricity 


f 


cost function 


fe 


fraction of the population doing local moves 


Id 


dominance index 


■^conv 


convergence metrics 




spreading metrics 




neighborhood of solution x 


Ne 


maximum number of allowed function evaluations 




number of function evaluations 




population size 


Pi 


i-th planet 


Pk 


population at generation k 


Pconv 


percentage of success on convergence 


Pspr 


percentage of success on spreading 


r 


random number 






'^pmin 


minimum pericentre radius 


s 


selection function 


s 


resource index 


T 


transfer time 


to 


departure time 


ti 


manoeuvre time 


tf 


final time 


U 


uniform distribution 
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u variation of the solution x 

X Pareto optimal set 

X solution vector 

y mutate individual 

Wc tolerance on the maximum allowable crowding 

Greek symbols 

Av variation of velocity 

p size of the neighborhood Np 

He gravity constant 

6i true anomaly of manoeuvre i 
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Table 1 : Multiobjective lest functions 



Scha 



f2 = {x- 5)2 







r 


if 


x<l 


X e [-5, 10] 


H 


-2 + x 

A- X 


^f 
^f 


1< a; < 3 

3 < a; < 4 






-4 + x 


if 


a; > 4 



Deb fi = xi 

a;i,a;2e[0,l] /2 = (l + 10a;2)[l-(^^)" 



a = 2; g = 4 


i+iox, sin(27rgxi) 




Deb2 


/l = Xl 




xi G [0, 1] 


f2 = g{Xi,X2)h(xi 


,x2);g{xi,x2) = 11 + a:2 — 10 cos(27ra;2) 


X2 e [-30, 30] 


h(xi,X2) = { ^ 


if /i < 9 


otherwise 


ZDT2 






Xi e [0,1]; 






i = 1, . . . ,n 


fi = a;i;/2 = gh 




n = 30 






ZDT4 


g = l + 10(n - 1) + E'Ui^i - 10cos(27rga;i)]; 


X, G [0,1]; 






Xi G [-5,5]; 


/i = Xi;f2 = gh 




i = 2, . . . ,n 






n = 10 






ZDT6 






Xi e [0,1]; 






i = 1, . . . ,n 


/i = 1 - exp(-4a;i) sin (Qnxi); f2 = gh 


n = 10 







Table 2: Comparison of the average Euclidean distances between 500 uniformly 
space points on the optimal Pareto front for various optimization algorithms: 
smaller dimension test problems. 



Approach 


Deb2 


Scha 


Deb 


MACS 


1.542e-3 


3.257e-3 


7.379e-4 




(5.19e-4) 


(5.61e-4) 


(6.36e-5) 


NSGA-II 


0.094644 


0.001594 


0.002536 




(0.117608) 


(0.000122) 


(0.000138) 


PAES 


0.259664 


0.070003 


0.002881 




(0.573286) 


(0.158081) 


(0.00213) 


MOPSO 


0.0011611 


0.00147396 


0.002057 




(0.0007205) 


(0.00020178) 


(0.000286) 
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Table 3: Comparison of the average Euclidean distances between 500 uniformly 
space points on the optimal Pareto front for various optimization algorithms: 
larger dimension test problems. 



Approach 


ZDT2 


ZDT4 


ZDT6 


MACS 


9.0896e-4 


0.0061 


0.0026 




(4.0862e-5) 


(0.0133) 


(0.0053) 


NSGA-II 


0.000824 


0.513053 


0.296564 




«le-5) 


(0.118460) 


(0.013135) 


PAES 


0.126276 


0.854816 


0.085469 




(0.036877) 


(0.527238) 


(0.006644) 



Table 4: Comparison of different versions of MACS. 



Approach 


Melric 


ZDT2 


ZDT4 


ZDT6 


Scha 


Deb 


Deb2 


MACS 


Pconv 


83.5% 


75% 


77% 


73% 


70.5% 


60% 




Pspr 


22.5% 


28% 


58.5% 


38.5% 


83% 


67.5 % 


MACS 


Pconv 


14% 


0% 


45% 


0.5% 


72.5% 


11% 


no local 


Pspr 


1% 


0% 


34% 


0% 


4% 


15% 


MACS 


Pconv 


84% 


22% 


78% 


37% 


92% 


21% 


p=l 


Pspr 


22% 


7% 


63% 


0% 


54% 


38% 


MACS 


Pconv 


56% 


42% 


57% 


78% 


85% 


42% 


p = 0.1 


Pspr 


5% 


15% 


61% 


88% 


94.5% 


74% 


MACS 


Pconv 


21% 


0.5% 


14.5% 


0.5% 


88.5% 


0% 


no attraction 


Pspr 


0% 


0.5% 


78.5% 


0% 


0% 


0% 



Table 5: Indexes Pconv and Pspr for different settings of MACS on the 3-impulse 
case 



pconv 






'^pop — 10 


^pop — 15 


fe = 


1/3 


45.5% 


55.5% 


61.0% 


fe = 


1/2 


48.0% 


51.0% 


55.5% 


fe = 


2/3 


45.0% 


52.5% 


43.0% 


Pspr 




TlpQp — 5 


^pop — 10 


^pop — 15 


/e = 


1/3 


68.5% 


62.0% 


56.0% 


fe = 


1/2 


65.0% 


57.0% 


46.0% 


fe = 


2/3 


67.5% 


51.0% 


36.5% 
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Table 6: NSGAII tuning on the 3-impulse case 



Mean Mconv 








Var Mconv 








Vc/Vm 


5 


25 


50 




5 


25 


50 


5 


36.1 


38.3 


43.0 


5 


201.0 


202.0 


185.0 


10 


32.3 


39.4 


40.6 


10 


182.0 


172.0 


182.0 


20 


31.7 


39.6 


42.5 


20 


175.0 


183.0 


169.0 


Mean M^pr 








Var Mgpr 










5 


25 


50 


Vc/Vm 


5 


25 


50 


5 


6.77 


7.24 


8.08 


5 


9.97 


9.47 


7.25 


10 


5.91 


7.50 


7.81 


10 


9.74 


8.68 


8.34 


20 


5.78 


7.50 


8.16 


20 


9.75 


8.53 


8.04 


Pconv 








Pspr 








Vc/Vm 


5 


25 


50 


•qj rim 


5 


25 


50 


5 


0.0% 


0.0% 


0.0% 


5 


44.8% 


37.7% 


23.4% 


10 


0.0% 


0.0% 


0.0% 


10 


57.8% 


33.1% 


29.2% 


20 


O.O'/r 


0.0% 


0.0%' 


20 


61.0% 


3 2. 5%' 


24.7% 



Table 7: MOPSO tuning on the 3-impulse case. 



Mean M^onv 








Var Mconv 








Particles/Subdivisions 


10 


30 


50 


Particles/Subdivisions 


10 


30 


50 


50 


59.1 


50.2 


41.5 


50 


1080.0 


1010.2 


713.3 


100 


50.0 


43.3 


41.3 


100 


591.0 


721.1 


778.1 


150 


47.8 


41.4 


39.4 


150 


562.1 


550.2 


608.2 


Mean Adgpr 








Var Mspr 








Particles/Subdivisions 


10 


30 


50 


Particles/Subdivisions 


10 


30 


50 


50 


16.1 


13.5 


12.3 


50 


41.3 


37.3 


24.0 


100 


14.7 


12.2 


11.6 


100 


32.8 


25.8 


24.5 


150 


14.8 


11.9 


11.4 


150 


30.0 


22.2 


22.5 


Pconv 








Pspr 








Particles/Subdivisions 


10 


30 


50 


Particles/Subdivisions 


10 


30 


50 


50 


0% 


0% 


0% 


50 


0.5% 


2.5% 


3.5% 


100 


0% 


0% 


0% 


100 


0.5% 


4.5% 


4% 


150 


0% 


0% 


0% 


150 


0.5% 


2% 


4% 
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Table 8: PAES tuning on the 3-impulse case 



Mean M^onv 








Var Mconv 








OUUUlVlMUllla/ iVlUullUJll 






n Q 


k> u uui V 1 aluii oi ivi u uiiiijii 




8 
w.o 


Q 


1 


53.7 


70.6 


70.2 


1 


525.0 


275.0 


297.0 


2 


52.8 


70.2 


70.0 


2 


479.0 


266.0 


305.0 


4 


53.0 


70.2 


70.1 


4 


453.0 


266.0 


311.0 


Mean Mgpr 








Var Mspr 








Subdivisions/Mutation 


0.6 


0.8 


0.9 


Subdivisions/Mutation 


0.6 


0.8 


0.9 


1 


14.2 


27.7 


36.7 


1 


20.0 


14.3 


17.3 


2 


13.6 


27.6 


36.6 


2 


17.5 


14.6 


16.8 


4 


13.8 


27.7 


36.6 


4 


17.2 


15.8 


17.0 


Pconv 
















Subdivisions/Mutation 


0.6 


0.8 


0.9 


Subdivisions/Mutation 


0.6 


0.8 


0.9 


1 


0.0% 


0.0% 


0.0% 


1 


0.0% 


0.0% 


0.0% 


2 


0.0% 


0.0% 


0.0% 


2 


0.0% 


0.0% 


0.0% 


4 


0.0% 


0.0% 


0.0% 


4 


0.0% 


0.0% 


0.0% 



Table 9: MTS tuning on the 3-impulse 



3imp 


Population 


20 


40 


80 


^conv 


Mean 


17.8 


22.6 


23.6 




Var 


97.6 


87.8 


73.2 




Mean 


12.6 


19.9 


18.4 




Var 


34.7 


26.2 


18.6 




Pconv 


1.0% 


0.0% 


0.0% 




Pspr 


0.5% 


0.0% 


0.0% 



Table 10: Summary of metrics Mconv and Mgpr, and associated performance 
indexes Pconv and Pspr, on the three impulse test cases. 



Metric 


MACS 


MACS 


MACS 


NSGA-II 


PAES 


MOPSO 


MTS 






no local 


no att 










■^conv 


5.53 


7.58 


154.7 


31.7 


53.0 


39.4 


17.8 




(15.1) 


(26.3) 


(235.0) 


(175.0) 


(453.0) 


(608.1) 


(97.6) 




5.25 


6.03 


9.16 


5.78 


13.8 


11.4 


12.6 




(3.73) 


(3.95) 


(2.07) 


(9.75) 


(17.2) 


(22.5) 


(34.7) 


Peon V 


61.0% 


40.5% 


0.0% 


0.0% 


0.0% 


0% 


1.0% 


Psp I- 


56.0% 


36Vr 


0.0% 


61.0% 


0.0% 


4.0% 


0.5% 
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Table 11: Metrics Mconv and M^pr, and associated performance indexes p, 
and Pspr, on the two impulse test cases. 



Metric MACS MACS MACS NSGA-II PAES MOPSO MTS 
no local no att 



■^conv 


0.0077 


0.0039 


0.0534 


0.283 


0.198 


0.378 


0.151 




(le-3) 


(1.4e-2) 


(9.5e-3) 


(4.35e-3) 


(0.332) 


(0.0636) 


(0.083) 




2.89 


6.87 


3.08 


2.47 


332.0 


2.11 


1.95 


(0.49) 


(7.36) 


0.943 


(0.119) 


(2.61e4) 


(2.65) 


(1.41) 


Pconv 


98.5% 


91.5% 


84% 


0% 


75.5% 


9% 


57.5% 


Pspr 


29.5% 


0% 


24% 


62.5% 


0.5% 


94.5% 


85.5% 
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Table 12: Metrics Mcmv and Mgpr, and associated peri'omiance indexes p, 
and Pspr, on the Cassini case. 



Approach 


Mctnc 


loUK 


jUUK 


OUUK 


MACS 




6.50 (229.1) 


4.48 (107.1) 


3.91 (62.6) 


5/15 




12.7 (126.0) 


11.1 (83.1) 


8.64 (35.5) 




Pconv 


6% 


14% 


21.5% 




Pspr 


Z/.jyo 


J 1 70 




MACS 




6.74 (217.9) 


5.56 (136.4) 


3.14 (50.1) 


5/10 




12.1 (83.0) 


10.3 (63.7) 


8.11 (30.0) 




Pconv 


8% 


10% 


25.5% 




Pspr 


ZO.Uyo 


Jl.j /c 


AC COI, 

43.370 


MACS 


■^conv 


13.5 (436.2) 


10.2 (350.1) 


7.86 (190.2) 


no local 




31.1 (274.3) 


27.9 (278.9) 


21.9 (226.7) 




Pconv 


1.0% 


1.5% 


2.5% 




Pspr 


1 A) /c 


L . J Vr 




MACS 


■^conv 


2.62 (1.46) 


2.2 (0.88) 


1.82 (0.478) 


no att 




27.8 (83.2) 


22.9 (58.0) 


18.2 (28.0) 




Pconv 


0.0% 


0.0% 


1.0% 




Pspr 






yj.yj /o 


NSGA-II 


■^conv 


2.43(18.0) 


1.99(16.8) 


1.24(1.62) 




Mgpr 


11.6 (71.4) 


11.0(47.5) 


8.78 (28.2) 




Pconv 


17.5% 


24.0% 


29.0% 




Pspr 


15.5% 


12.5% 


25.0% 


MOPSO 


■^conv 


2.62 (7.33) 


2.4 (2.57) 


2.14(0.94) 






28.0 (308.3) 


24.6 (260.4) 


21.8 (231.3) 




Pconv 


0.5% 


1.0% 


1.0% 




Pspr 


0.0% 


0.0% 


0.5% 


PAES 


■^conv 


24.0 (54.5) 


19.8 (32.9) 


15.2(16.6) 




]\d^gpf 


30.1 (47.3) 


26.0 (33.9) 


21.4(19.5) 




Pconv 


0.0% 


0.0% 


0.0% 




Pspr 


0.0% 


0.0% 


0.0% 


MTS 


^conv 


3.71 (1.53) 


3.39(1.67) 


3.02(1.69) 




^spr 


18.1 (18.2) 


15.6(13.4) 


13.1 (8.46) 




Pconv 


0.0% 


0.0% 


0.0% 




Pspr 


0.0% 


0.0% 


0.0% 
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Figure 2: Three-impulse test case: a) Complete Pareto front, b) close-up of the 
Pareto fronts 
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Figure 3: Cassini test case: a) Complete Pareto front, b) close-up of the Pareto 
fronts 
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